#------------------------------------------------------------------------------
# include all functions to generate Figure 1,2,3 and S3 in SI.
#==============================================================================

trim <- function (x) gsub("^\\s+|\\s+$", "", x)

list_target_category = list()
list_target_category[[1]] = list(name='Female',value=c(0,1),value_name=c('Male','Female'))
list_target_category[[2]] = list(name='AgeGrp4',value=c(1,2,3,4),value_name=c('15_24','25_44','45_64','65_Up'))
list_target_category[[3]] = list(name='Race5',value=c(1,2,3,4,5),value_name=c('White','Black','Indian|Alaska','Asian','Hispanic'))
list_target_category[[4]] = list(name='MarStat5',value=c(1,2,3,4,5),value_name=c('Married','Widowed','Divorced','Separated','Never Married'))
list_target_category[[5]] = list(name='BornUSA',value=c(0,1),value_name=c('not US-born','US-born'))
list_target_category[[6]] = list(name='UnEmpl',value=c(0,1),value_name=c('Employed','Unemployed'))
list_target_category[[7]] = list(name='UnEmpl_raw',value=c(0,1),value_name=c('Employed','Unemployed'))
list_target_category[[8]] = list(name='UnEmpl_emp',value=c(0,1),value_name=c('Employed','Unemployed'))
list_target_category[[9]] = list(name='UnEmpl_unemp',value=c(0,1),value_name=c('Employed','Unemployed'))
list_target_category[[10]] = list(name='UnEmpl_mean',value=c(0,1),value_name=c('Employed','Unemployed'))
list_target_category[[11]] = list(name='PhysProb',value=c(0,1),value_name=c('Healthy','Physical Problem'))
list_target_category[[12]] = list(name='PhysProb n~b',value=c(0,1),value_name=c('Healthy','Physical Problem'))
list_target_category[[13]] = list(name='PhysProb p~b',value=c(0,1),value_name=c('Healthy','Physical Problem'))
list_target_category[[14]] = list(name='PhysProb m~n',value=c(0,1),value_name=c('Healthy','Physical Problem'))

find_category = function(target){
	i = which(unlist(lapply(list_target_category,function(x) x$name)) %in% target)
	value = list_target_category[[i]]$value 
	value_name = list_target_category[[i]]$value_name
	return(list(value=value,value_name=value_name))
}

match_category = function(target){
	grepl(paste0(unlist(lapply(list_target_category,function(x) x$name)),collapse='|'),
			target)
}

match_continuous = function(target){
	grepl('Pop|Rat|RAT|prop|deviant|same',target)
}

read_output = function(log_path,method='logit',sameness='raw', mm=1){
	
	results = readLines(file.path(log_path,paste0(method,'_',mm,'.txt')))	
	
	results = trim(results)		

	# some targets checking
	split2 = c('-------------+-------------------------------------------------------')

	target_at = which(grepl('Changes in Pr(y)',results,fixed=TRUE))
	target_beg = which(results %in% split2)
	target_end = which(results %in% 'Average predictions')
	
	target_loc = NULL
	for (i in target_at) {
		beg = min(target_beg[target_beg>i])
		end = min(target_end[target_end>i])
		target_loc = rbind(target_loc,c(beg+1,end-1))
	}
	
	target_loc = target_loc[target_loc[,1] != Inf|target_loc[,2] != Inf,]
	
	target_output = results[target_loc[1]:target_loc[2]]	

	list_output1 = list(); n1 = 1
	list_output2 = list(); n2 = 1

	for (j in 1:length(target_output)){
		X = target_output[j]
		#message(j,X)
		if (match_category(X) & !match_continuous(X)) {
			category = trim(unlist(strsplit(X,'\\|')))
			category = category[category != '']
		}

		if (match_continuous(X)) {
			#if (grepl('deviant',X)) category = paste0('majority ',category)
			if (grepl('same',X)) {
				category = paste0('sameness ',category)

			} else {
			category = trim(unlist(strsplit(X,'\\|')))
			category = category[category != '']
			}
		}

		if (grepl('vs|to|\\+|delta|Range|Marginal',X)) {
			stat = unlist(strsplit(strsplit(X,'\\|')[[1]][2],' '))
			stat = stat[stat!='']
			if (length(stat) == 5) {
				names(stat) = c('change','p','ll','ul','z')	
				vs = trim(strsplit(X,'\\|')[[1]][1])
				output = data.table(category,mtype=vs,t(stat))
				list_output1[[n1]] = output; n1 = n1 + 1

			} else if (length(stat) == 3) {
				names(stat) = c('se','from','to')	
				vs = trim(strsplit(X,'\\|')[[1]][1])
				output = data.table(category,mtype=vs,t(stat))
				list_output2[[n2]] = output; n2 = n2 + 1
			}
		}
	}


	list_output1 = rbindlist(list_output1)	
	list_output2 = rbindlist(list_output2)	
	list_output = merge(x=list_output1, y=list_output2, by=c('category','mtype'),all.x=TRUE,all.y=TRUE)

	list_output[category == 'RAT Female', category := 'prop_Female']
	list_output[category == 'RAT AgeGrp~2', category := 'prop_AgeGrp4_2']
	list_output[category == 'RAT AgeGrp~3', category := 'prop_AgeGrp4_3']
	list_output[category == 'RAT AgeGrp~4', category := 'prop_AgeGrp4_4']
	list_output[category == 'RAT Race5 2', category := 'prop_Race5_2']
	list_output[category == 'RAT Race5 3', category := 'prop_Race5_3']
	list_output[category == 'RAT Race5 4', category := 'prop_Race5_4']
	list_output[category == 'RAT Race5 5', category := 'prop_Race5_5']
	list_output[category == 'RAT BornUSA', category := 'prop_BornUSA']	
	list_output[category == 'RAT MarSta~2', category := 'prop_MarStat5_2']
	list_output[category == 'RAT MarSta~3', category := 'prop_MarStat5_3']
	list_output[category == 'RAT MarSta~4', category := 'prop_MarStat5_4']
	list_output[category == 'RAT MarSta~5', category := 'prop_MarStat5_5']

	list_output[category == 'RAT UnEmpl', category := 'prop_UnEmpl']
	list_output[category == 'RAT UnEmpl~w', category := 'prop_UnEmpl']
	list_output[category == 'RAT UnEmpl~p', category := 'prop_UnEmpl']
	list_output[category == 'RAT UnEmpl~n', category := 'prop_UnEmpl']

	list_output[category == 'UnEmpl raw', category := 'UnEmpl']
	list_output[category == 'UnEmpl unemp', category := 'UnEmpl']
	list_output[category == 'UnEmpl emp', category := 'UnEmpl']
	list_output[category == 'UnEmpl mean', category := 'UnEmpl']

	list_output[category == 'sameness UnEmpl unemp', category := 'sameness UnEmpl']
	list_output[category == 'sameness UnEmpl emp', category := 'sameness UnEmpl']
	list_output[category == 'sameness UnEmpl raw', category := 'sameness UnEmpl']
	list_output[category == 'sameness UnEmpl mean', category := 'sameness UnEmpl']

	list_output[grep('RAT PhysPr',category), category := 'prop_PhysProb']
	list_output[category == 'RAT PhysProb', category := 'prop_PhysProb']
	list_output[category == 'PhysProb raw', category := 'PhysProb']
	list_output[category == 'PhysProb n~b', category := 'PhysProb']
	list_output[category == 'PhysProb p~b', category := 'PhysProb']
	list_output[category == 'PhysProb m~n', category := 'PhysProb']

	list_output[category == 'sameness PhysProb raw', category := 'sameness PhysProb']
	list_output[category == 'sameness PhysProb n~b', category := 'sameness PhysProb']
	list_output[category == 'sameness PhysProb p~b', category := 'sameness PhysProb']
	list_output[category == 'sameness PhysProb m~n', category := 'sameness PhysProb']

	coef_combine = data.table(method=method, model_type = mm, sameness = sameness, list_output)

	coef_combine$at_name = factor(coef_combine$category,
	levels=c(c('Female','AgeGrp4','Race5','BornUSA','MarStat5','UnEmpl','PhysProb'),
		'prop_Female','prop_AgeGrp4_2','prop_AgeGrp4_3','prop_AgeGrp4_4',
		'prop_Race5_2','prop_Race5_3','prop_Race5_4','prop_Race5_5',
		'prop_BornUSA','prop_MarStat5_2','prop_MarStat5_3','prop_MarStat5_4','prop_MarStat5_5','prop_UnEmpl',
		'Rat Poverty','Rat Mig Cum','Pop Den',
		'Rat GC ProE','Rat GC ProM','Rat GC ProB','Rat GC Cat~x','Rat GC Jew','RAT GC Mor',
		paste0('sameness ',c('Female','AgeGrp4','Race5','BornUSA','MarStat5','UnEmpl','PhysProb')),
		paste0('majority ',c('Female','AgeGrp4','Race5','BornUSA','MarStat5','UnEmpl','PhysProb'))
		),
	labels=c(
		c('Female','Age Group','Race','BornUSA','Marital Status','UnEmployed','Physical Problem'),
		paste0(c('% Female','% 25-44 age group','% 45-64 age group', '% 65+ age group',
		'% Black','% Indian|Alaska Native','% Asian','% Hispanic',
		'% US Born', '% Widowed','% Divorced','% Separated','% Never Married','% Unemployed'),' in PUMA Group'),
		'% Below Poverty Line','% Net Migration','Population Density',
		'Religion: % Evangelical protestant', 'Religion: % Mainline protestant', 'Religion: % Black protestant','Religion: % Catholic Orthdox','Religion: % Jewish','Religion: %Mormon',
		paste0('% Same ',c('Sex','Age Group','Race','BornUSA','Marital Status','UnEmployed','Physical Problem')),
		paste0('% Majority ',c('Sex','Age Group','Race','BornUSA','Marital Status','UnEmployed','Physical Problem'))))

	coef_combine$at_type =  gsub('Rat |prop_|sameness ','',coef_combine$category)
	coef_combine[grepl('AgeGrp4',at_type), at_type := 'Age']
	coef_combine[grepl('GC',at_type), at_type := 'Religion']
	coef_combine[grepl('MarStat5',at_type), at_type := 'Marriage']
	coef_combine[grepl('Race5',at_type), at_type := 'Race']
	coef_combine[grepl('Den|Poverty|Mig',at_type), at_type := 'Others']
	coef_combine[,at_type := factor(at_type, levels=c('Female','Age','Race','BornUSA',
		'Marriage','UnEmpl','PhysProb','Religion','Others'),
	labels=c('Sex','Age','Race','USA Born','Marital Status','Unemployment','Physical Problem','Religious Congregation','Other factors'))]

	coef_combine[,category_name := mtype]

	coef_combine[category == 'Female' & mtype == '1 vs 0',category_name := 'Female (vs Male)']
	coef_combine[category == 'prop_Female' & mtype == 'Range',category_name := '% Female [range]']

	coef_combine[category == 'AgeGrp4' & mtype == '4 vs 1',category_name := 'Age 65+ (vs Age 15-24)']
	coef_combine[category == 'AgeGrp4' & mtype == '3 vs 1',category_name := 'Age 45-64 (vs Age 15-24)']
	coef_combine[category == 'AgeGrp4' & mtype == '2 vs 1',category_name := 'Age 25-44 (vs Age 15-24)']

	coef_combine[category == 'prop_AgeGrp4_2',category_name := '% Age 65+ [range]']
	coef_combine[category == 'prop_AgeGrp4_3',category_name := '% Age 45-64 [range]']
	coef_combine[category == 'prop_AgeGrp4_4',category_name := '% Age 25-44 [range]']


	coef_combine[category == 'Race5' & mtype == '2 vs 1',category_name := 'Black (vs White)']
	coef_combine[category == 'Race5' & mtype == '3 vs 1',category_name := 'American Indian (vs White)']
	coef_combine[category == 'Race5' & mtype == '4 vs 1',category_name := 'Asian (vs White)']
	coef_combine[category == 'Race5' & mtype == '5 vs 1',category_name := 'Hispanic (vs White)']

	coef_combine[category == 'prop_Race5_2',category_name := '% Black [range]']
	coef_combine[category == 'prop_Race5_3',category_name := '% American Indian [range]']
	coef_combine[category == 'prop_Race5_4',category_name := '% Asian [range]']
	coef_combine[category == 'prop_Race5_5',category_name := '% Hispanic [range]']

	coef_combine[category == 'BornUSA' & mtype == '1 vs 0',category_name := 'Native (vs non-native)']
	coef_combine[category == 'prop_BornUSA' & mtype == 'Range',category_name := '% Native [range]']

	coef_combine[category == 'MarStat5' & mtype == '2 vs 1',category_name := 'Widowed (vs Married)']
	coef_combine[category == 'MarStat5' & mtype == '3 vs 1',category_name := 'Divorced (vs Married)']
	coef_combine[category == 'MarStat5' & mtype == '4 vs 1',category_name := 'Separated (vs Married)']
	coef_combine[category == 'MarStat5' & mtype == '5 vs 1',category_name := 'Never-married (vs Married)']

	coef_combine[category == 'prop_MarStat5_2',category_name := '% Widowed [range]']
	coef_combine[category == 'prop_MarStat5_3',category_name := '% Divorced [range]']
	coef_combine[category == 'prop_MarStat5_4',category_name := '% Separated [range]']
	coef_combine[category == 'prop_MarStat5_5',category_name := '% Never-married [range]']
	
	coef_combine[category == 'UnEmpl' & mtype == '1 vs 0',category_name := 'Unemployed (vs employed)']
	coef_combine[category == 'UnEmpl' & mtype == 'Range',category_name := 'Unemployed (vs employed)']

	coef_combine[category == 'prop_UnEmpl' & mtype == 'Range',category_name := '% Unemployed [range]']

	coef_combine[category == 'PhysProb' & mtype == '1 vs 0',category_name := 'Physical Problem (vs No Problem)']
	coef_combine[category == 'PhysProb' & mtype == 'Range',category_name := 'Physical Problem (vs No Problem)']

	coef_combine[category == 'prop_PhysProb' & mtype == 'Range',category_name := '% Physical Problem [range]']

	coef_combine[category == 'Rat GC ProE',category_name := '% Evangelical protestant [range]']
	coef_combine[category == 'Rat GC ProM',category_name := '% Mainline protestant [range]']
	coef_combine[category == 'Rat GC ProB' ,category_name := '% Black protestant [range]']
	coef_combine[category == 'Rat GC Cat~x',category_name := '% Catholic Orthdox [range]']
	coef_combine[category == 'Rat GC Jew'  ,category_name := '% Jewish [range]']

	coef_combine[category == 'Pop Den' & mtype == 'Range',category_name := '% population density [range]']
	coef_combine[category == 'Rat Poverty' & mtype == 'Range',category_name := '% below poverty [range]']
	coef_combine[category == 'Rat Mig Cum' & mtype == 'Range',category_name := '% net migration [range]']

	coef_combine[sameness == 1 & mtype == 'Range' &!(at_type %in% c('Religious Congregation','Other factors')), category_name := '% Same Category [range]']

	coef_combine[,change := as.numeric(change)]
	coef_combine[,abs_change := abs(change)]
	coef_combine[,ll := as.numeric(ll)]
	coef_combine[,ul := as.numeric(ul)]

	return(coef_combine)
}


read_output_mi = function(log_path,method='mi',sameness='raw', mm=1){
	
	results = readLines(file.path(log_path,paste0(method,'_',mm,'.txt')))	
	
	results = trim(results)		
	split2 = c('-------------+------------------------------------------------------------')
	splitg = c('-------------+------------------------------------------------')

	if (mm != 9){
		target2 = which(results %in% split2) + 1
		output2 = results[target2]
	} else {
		target2 = c(which(results %in% split2) + 1,
		448:453, 478:487, 561:570, 641,666)
		output2 = results[target2]
		output2 = grep(' vs 0',output2,value=TRUE)
	}
	
	targetg_beg = which(results %in% c(splitg))

	if (mm == 9){
		names_output2 = c('UnEmpl','PhysProb')
	}  else if (mm == 10){
		names_output2 = c('prop_Female','prop_AgeGrp4_2','prop_AgeGrp4_3','prop_AgeGrp4_4',
			'prop_Race5_2','prop_Race5_3','prop_Race5_4','prop_Race5_5','prop_BornUSA',
			'prop_MarStat5_2','prop_MarStat5_3','prop_MarStat5_4','prop_MarStat5_5','prop_UnEmpl','prop_PhysProb')
	} else {
	if (sameness == 'raw'){
		names_output2 = c('Female','BornUSA','prop_Female','prop_AgeGrp4_2','prop_AgeGrp4_3','prop_AgeGrp4_4',
			'prop_Race5_2','prop_Race5_3','prop_Race5_4','prop_Race5_5','prop_BornUSA',
			'prop_MarStat5_2','prop_MarStat5_3','prop_MarStat5_4','prop_MarStat5_5','prop_UnEmpl','prop_PhysProb')
	} else if (sameness == 'sameness'){
		names_output2 = c('Female','BornUSA',
			paste0('sameness ',c('Female','AgeGrp4','Race5','BornUSA','MarStat5','UnEmpl','PhysProb')))
	}		
	}

	list_output2 = list()
	for (i in 1:length(output2)){
		tmp = unlist(strsplit(strsplit(output2[i],'\\| ')[[1]][2],' '))
		tmp = tmp[tmp!='']
		if (mm != 9){
			names(tmp) = c('coef','se','z','p','ll','ul')	
			list_output2[[i]] = data.frame(t(tmp[c(1,2,5,6)]), category = names_output2[i])
		} else {
			names(tmp) = c('coef','se','ll','ul')	
			list_output2[[i]] = data.frame(t(tmp), category = names_output2[i])			
		}		
	}
	list_output2 = rbindlist(list_output2)

	if (mm == 9) {
		list_output2 = data.frame(category_name = rep('1 vs 0',2),list_output2)
	} else if (mm == 10){
		list_output2 = data.frame(category_name = rep('Range',15),list_output2)
	} else {
		list_output2 = data.frame(category_name = c('1 vs 0','1 vs 0',rep('Range',nrow(list_output2)-2)),list_output2)
	}

	if (mm == 9 | mm == 10){

		list_output = list_output2

	} else {

	list_output3 = list(); k = 1
	for (i in 1:length(targetg_beg)){

		beg = targetg_beg[i] + 2 
		if (i == 1) {
			end = beg + 2 
			category = 'AgeGrp4'
		}
		if (i == 2) {
			end = beg + 3 
			category = 'Race5'			
		}
		if (i == 3) {
			end = beg + 3 
			category = 'MarStat5'
		}

		if (i == 4) {
			end = beg 
			category = 'UnEmpl'
		}

		if (i == 5) {
			end = beg 
			category = 'PhysProb'
		}

		tmp_all = results[beg:end]
		list_tmp = list()
		for (j in 1:length(tmp_all)){
			tmp = tmp_all[j]
			name = trim(strsplit(tmp,'\\| ')[[1]][1])
			tmp = unlist(strsplit(strsplit(tmp,'\\| ')[[1]][2],' '))
			tmp = tmp[tmp!='']
			tmp = c(name,tmp)
			names(tmp) = c('category_name','coef','se','ll','ul')
			list_tmp[[j]] = data.frame(t(tmp))
		}
		list_tmp = rbindlist(list_tmp)
		list_output3[[k]] = data.frame(list_tmp,category=category)
		k = k + 1
	}
	list_output3 = rbindlist(list_output3)
	
	list_output = rbind(list_output2,list_output3)		
	}

	coef_combine = data.table(method=method, model_type = mm, sameness = sameness, list_output)

	coef_combine$at_name = factor(coef_combine$category,
	levels=c(c('Female','AgeGrp4','Race5','BornUSA','MarStat5','UnEmpl','PhysProb'),
		'prop_Female','prop_AgeGrp4_2','prop_AgeGrp4_3','prop_AgeGrp4_4',
		'prop_Race5_2','prop_Race5_3','prop_Race5_4','prop_Race5_5',
		'prop_BornUSA','prop_MarStat5_2','prop_MarStat5_3','prop_MarStat5_4','prop_MarStat5_5','prop_UnEmpl',
		'Rat Poverty','Rat Mig Cum','Pop Den',
		'Rat GC ProE','Rat GC ProM','Rat GC ProB','Rat GC Cat~x','Rat GC Jew','RAT GC Mor',
		paste0('sameness ',c('Female','AgeGrp4','Race5','BornUSA','MarStat5','UnEmpl','PhysProb')),
		paste0('majority ',c('Female','AgeGrp4','Race5','BornUSA','MarStat5','UnEmpl','PhysProb'))
		),
	labels=c(
		c('Female','Age Group','Race','BornUSA','Marital Status','UnEmployed','Physical Problem'),
		paste0(c('% Female','% 25-44 age group','% 45-64 age group', '% 65+ age group',
		'% Black','% Indian|Alaska Native','% Asian','% Hispanic',
		'% US Born', '% Widowed','% Divorced','% Separated','% Never Married','% Unemployed'),' in PUMA Group'),
		'% Below Poverty Line','% Net Migration','Population Density',
		'Religion: % Evangelical protestant', 'Religion: % Mainline protestant', 'Religion: % Black protestant','Religion: % Catholic Orthdox','Religion: % Jewish','Religion: %Mormon',
		paste0('% Same ',c('Sex','Age Group','Race','BornUSA','Marital Status','UnEmployed','Physical Problem')),
		paste0('% Majority ',c('Sex','Age Group','Race','BornUSA','Marital Status','UnEmployed','Physical Problem'))))

	coef_combine$at_type =  gsub('Rat |prop_|sameness ','',coef_combine$category)
	coef_combine[grepl('AgeGrp4',at_type), at_type := 'Age']
	coef_combine[grepl('GC',at_type), at_type := 'Religion']
	coef_combine[grepl('MarStat5',at_type), at_type := 'Marriage']
	coef_combine[grepl('Race5',at_type), at_type := 'Race']
	coef_combine[grepl('Den|Poverty|Mig',at_type), at_type := 'Others']
	coef_combine[,at_type := factor(at_type, levels=c('Female','Age','Race','BornUSA',
		'Marriage','UnEmpl','PhysProb','Religion','Others'),
	labels=c('Sex','Age','Race','USA Born','Marital Status','Unemployment','Physical Problem','Religious Congregation','Other factors'))]

	setnames(coef_combine,'category_name','mtype')
	coef_combine[,category_name := mtype]

	coef_combine[category == 'Female' & mtype == '1 vs 0',category_name := 'Female (vs Male)']
	coef_combine[category == 'prop_Female' & mtype == 'Range',category_name := '% Female [range]']

	coef_combine[category == 'AgeGrp4' & mtype == '4 vs 1',category_name := 'Age 65+ (vs Age 15-24)']
	coef_combine[category == 'AgeGrp4' & mtype == '3 vs 1',category_name := 'Age 45-64 (vs Age 15-24)']
	coef_combine[category == 'AgeGrp4' & mtype == '2 vs 1',category_name := 'Age 25-44 (vs Age 15-24)']

	coef_combine[category == 'prop_AgeGrp4_2',category_name := '% Age 65+ [range]']
	coef_combine[category == 'prop_AgeGrp4_3',category_name := '% Age 45-64 [range]']
	coef_combine[category == 'prop_AgeGrp4_4',category_name := '% Age 25-44 [range]']


	coef_combine[category == 'Race5' & mtype == '2 vs 1',category_name := 'Black (vs White)']
	coef_combine[category == 'Race5' & mtype == '3 vs 1',category_name := 'American Indian (vs White)']
	coef_combine[category == 'Race5' & mtype == '4 vs 1',category_name := 'Asian (vs White)']
	coef_combine[category == 'Race5' & mtype == '5 vs 1',category_name := 'Hispanic (vs White)']

	coef_combine[category == 'prop_Race5_2',category_name := '% Black [range]']
	coef_combine[category == 'prop_Race5_3',category_name := '% American Indian [range]']
	coef_combine[category == 'prop_Race5_4',category_name := '% Asian [range]']
	coef_combine[category == 'prop_Race5_5',category_name := '% Hispanic [range]']

	coef_combine[category == 'BornUSA' & mtype == '1 vs 0',category_name := 'Native (vs non-native)']
	coef_combine[category == 'prop_BornUSA' & mtype == 'Range',category_name := '% Native [range]']

	coef_combine[category == 'MarStat5' & mtype == '2 vs 1',category_name := 'Widowed (vs Married)']
	coef_combine[category == 'MarStat5' & mtype == '3 vs 1',category_name := 'Divorced (vs Married)']
	coef_combine[category == 'MarStat5' & mtype == '4 vs 1',category_name := 'Separated (vs Married)']
	coef_combine[category == 'MarStat5' & mtype == '5 vs 1',category_name := 'Never-married (vs Married)']

	coef_combine[category == 'prop_MarStat5_2',category_name := '% Widowed [range]']
	coef_combine[category == 'prop_MarStat5_3',category_name := '% Divorced [range]']
	coef_combine[category == 'prop_MarStat5_4',category_name := '% Separated [range]']
	coef_combine[category == 'prop_MarStat5_5',category_name := '% Never-married [range]']
	
	coef_combine[category == 'UnEmpl' & mtype == '1 vs 0',category_name := 'Unemployed (vs employed)']
	coef_combine[category == 'UnEmpl' & mtype == 'Range',category_name := 'Unemployed (vs employed)']

	coef_combine[category == 'prop_UnEmpl' & mtype == 'Range',category_name := '% Unemployed [range]']

	coef_combine[category == 'PhysProb' & mtype == '1 vs 0',category_name := 'Physical Problem (vs No Problem)']
	coef_combine[category == 'PhysProb' & mtype == 'Range',category_name := 'Physical Problem (vs No Problem)']

	coef_combine[category == 'prop_PhysProb' & mtype == 'Range',category_name := '% Physical Problem [range]']

	coef_combine[category == 'Rat GC ProE',category_name := '% Evangelical protestant [range]']
	coef_combine[category == 'Rat GC Mor' ,category_name := '% Mormon [range]']
	coef_combine[category == 'Rat GC ProM',category_name := '% Mainline protestant [range]']
	coef_combine[category == 'Rat GC ProB' ,category_name := '% Black protestant [range]']
	coef_combine[category == 'Rat GC Cat~x',category_name := '% Catholic Orthdox [range]']
	coef_combine[category == 'Rat GC Jew'  ,category_name := '% Jewish [range]']

	coef_combine[category == 'Pop Den' & mtype == 'Range',category_name := '% population density [range]']
	coef_combine[category == 'Rat Poverty' & mtype == 'Range',category_name := '% below poverty [range]']
	coef_combine[category == 'Rat Mig Cum' & mtype == 'Range',category_name := '% net migration [range]']

	coef_combine[sameness == 1 & mtype == 'Range' &!(at_type %in% c('Religious Congregation','Other factors')), category_name := '% Same Category [range]']

	coef_combine[,change := as.numeric(as.character(coef))]
	coef_combine[,ll := as.numeric(as.character(ll))]
	coef_combine[,ul := as.numeric(as.character(ul))]

	return(coef_combine)
}

read_interaction_plot = function(outcome, target_var, 
	method='logit',sameness='sameness', model_no=5,log_path=log_path){

	results = readLines(file.path(log_path,paste0(method,'_',model_no,'.txt')))	
	results = trim(results)

# some targets checking
split2 = c(
	'--------------------------------------------------------------------------------',
	'------------------------------------------------------------------------------',
	'----------------------------------------------------------------------------------')

target_at = which(grepl('_at#',results))
target_end = which(results %in% split2)

target_loc = NULL
for (i in target_at) {
	j = min(target_end[target_end>i])
	target_loc = rbind(target_loc,c(i+1,j-1))
}

target_at = which(results %in% 'Expression   : Pr(Suic), predict(pr)')
target_end = which(results %in% split2)

target_loc_at = NULL
for (i in target_at) {
	j = min(target_end[target_end>i])
	target_loc_at = rbind(target_loc_at,c(i+1,j-1))
}

target_names = trim(gsub('*_at\\#| \\|','',grep('_at#',results,value=TRUE)))

list_coef_tab = list()

for (i in 1:length(target_names)){
	target = target_names[i]
	target_by = find_category(target)

	target_at = results[target_loc_at[i,1]:target_loc_at[i,2]]
	target_at = strsplit(gsub('\\._at|\\:|\\=','',target_at),' ')
	target_at = lapply(target_at,function(x) x[x!=''])

	target_at = as.data.frame(do.call(rbind,target_at))
	names(target_at) = c('at','varname','value')
	target_at$value = as.numeric(as.character(target_at$value))
	target_at$at = as.numeric(as.character(target_at$at))

	coef_tab = results[target_loc[i,1]:target_loc[i,2]]
	coef_tab = strsplit(gsub('\\|',' ', coef_tab), ' ')
	
	coef_tab_out = list()
	for (j in 1:length(coef_tab)) {
		coef = coef_tab[[j]]
		coef = coef[coef != '']
		if (length(coef) == 8) {
			coef_dt = data.table(t(as.numeric(coef)))
		} 
		if (length(coef) == 5) {
			coef_dt = data.table(t(c(na.omit(as.numeric(coef)),rep(NA,6))))
		}
		coef_tab_out[[j]] = coef_dt
	}

	coef_tab_out = rbindlist(coef_tab_out)
	colnames(coef_tab_out) = c('at','ind','margin','se','t','p','lci','uci')

	coef_tab_out$at_value = target_at$value[match(coef_tab_out$at,target_at$at)]
	coef_tab_out$at_name = target_at$varname[match(coef_tab_out$at,target_at$at)]

	coef_tab_out$ind_value = factor(target_by$value_name[match(coef_tab_out$ind,target_by$value)],
		levels=target_by$value_name)
	coef_tab_out$var_name = target 
	coef_tab_out$model_no = model_no
	list_coef_tab[[i]] = coef_tab_out
}	

coef_combine = rbindlist(list_coef_tab)
return(coef_combine)
}

plot_interaction = function(outcome, target_var, method='logit',sameness='sameness', model_no=3, 
	log_path=log_path,plot_scale ='free_x'){

	coef_combine = read_interaction_plot(outcome=outcome,target_var=target_Var,
		method=method,sameness=sameness,model_no=model_no,log_path=log_path)
	
	coef_combine$at_name = factor(coef_combine$at_name,
		levels=c('std_same_p~x','std_same_p~4','std_same_~e5','std_same_p~A','std_same_~t5','std_same_p~l','std_same_p~b'),
		labels=c(paste0('Standardized ',c('% Same Sex','% Same Age Group','% Same Race','% Same US Born', '% Same Martial Status', '% Same Employment Status','% Same Physical Problem'),'')))
	
	coef_for_plot = coef_combine %>% filter(var_name %in% target_var & model_no == model_no)
	
	legend_name = gsub('\\d','',target_var)
	
	if (legend_name == 'Female') legend_name <- 'Sex'
	if (legend_name == 'MarStat') legend_name <- 'Marital Status'
	if (grepl('UnEmpl',legend_name)) legend_name <- 'Employment Status'
	if (legend_name == 'BornUSA') legend_name <- 'USA Born'
	if (legend_name == 'AgeGrp') legend_name <- 'Age Group'
	
	p <- coef_for_plot %>% 
		ggplot(aes(x=at_value,y=margin*100000,ymin=lci*100000,ymax=uci*100000,
		shape=ind_value,color=ind_value,fill=ind_value)) +
	 	geom_ribbon(alpha=0.3,size=0)+
	 	#geom_pointrange(lty=1, position=position_dodge(width=1))+
	 	#geom_pointrange(lty=1,lwd=0.5)+
	 	geom_point(size=2)+
	 	geom_line() + theme_pubr() +
	 	scale_fill_nejm()+
	 	scale_color_nejm()+
	 	theme(legend.position = "top",
	 		legend.text=element_text(size=5),
	 		legend.key.width=unit(0.3, "cm"),
	 		legend.key.height=unit(0.3, "cm"),
	 		legend.title=element_blank()) +
		labs(title = NULL, y="Suicide rate per 100K population")
		
	if (any(legend_name %in% c('Marital Status','Race'))){
		p = p + guides(
			color=guide_legend(title=legend_name, nrow=2),
			fill=guide_legend(title=legend_name,nrow=2),
			shape=guide_legend(title=legend_name,nrow=2)) 
	} else {
		p = p + guides(
			color=guide_legend(title=legend_name),
			fill=guide_legend(title=legend_name),
			shape=guide_legend(title=legend_name)) 
	}
	
	p <- p + labs(x=coef_for_plot$at_name, subtitle=NULL)
	
	
	return(p)
}



